Testing the different components of the enveloppe
Motion Clouds: testing components of the envelope¶
In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
import os, sys
sys.path.append('..')
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
#mc.recompute = True
mc.notebook = True
Testing the color¶
In [3]:
%bookmark root
In [4]:
%cd..
In [5]:
%run test_color.py
Testing the speed¶
In [6]:
%run test_speed.py
Exploring the orientation component of the envelope around a grating.¶
In [7]:
%run test_grating.py
In [8]:
%run test_radial.py
In [9]:
%cd -b root
Testing competing Motion Clouds with psychopy
Analysis of a discrimination experiment¶
In the psychopy_competition.py script, we implemented an experiment to test whether one could discriminate the dominating motionin a sum of two motion cloudsin opposite directions.
Herein, we analyse the data that was collected over different sessions and try to draw some conclusions.
In []:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
The experiment¶
In [1]:
#%pycat psychopy_competing.py
Analysis, version 1¶
In a first version of the experiment, we only stored the results in a numpy file:
In [2]:
!ls data/competing_v1_*npy
In [3]:
import glob
fig = plt.figure()
ax = fig.add_subplot(111)
for fn in glob.glob('data/competing_v1_*npy'):
results = np.load(fn)
print 'experiment ', fn, ', # trials=', results.shape[1]
ax.scatter(results[1, :], results[0, :])
_ = ax.axis([0., 1., -1.1, 1.1])
_ = ax.set_xlabel('contrast')
_ = ax.set_ylabel('perceived direction')
In [4]:
alpha = .3
fig = plt.figure(figsize=(12,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
data = []
for fn in glob.glob('data/competing_v1_*npy'):
results = np.load(fn)
data_ = np.empty(results.shape)
# lower stronger, detected lower = CORRECT is 1
ax1.scatter(results[1, results[1, :]<.5],
1.*(results[0, results[1, :]<.5]==-1),
c='r', alpha=alpha)
# copying data: contrast (from .5 to 1), correct
data_[0, results[1, :]<.5] = 1. - results[1, results[1, :]<.5]
data_[1, results[1, :]<.5] = 1.*(results[0, results[1, :]<.5]==-1)
# upper stronger, detected lower = INCORRECT is 1
ax1.scatter(results[1, results[1, :]>.5],
1.*(results[0, results[1, :]>.5]==-1),
c='b', alpha=alpha)
# lower stronger, detected upper = INCORRECT is 1
ax2.scatter(results[1, results[1, :]<.5],
1.*(results[0, results[1, :]<.5]==1),
c='b', alpha=alpha, marker='x')
# upper stronger, detected upper = CORRECT is 1
ax2.scatter(results[1, results[1, :]>.5],
1.*(results[0, results[1, :]>.5]==1),
c='r', alpha=alpha, marker='x')
# copying data: contrast (from .5 to 1), correct
data_[0, results[1, :]>=.5] = results[1, results[1, :]>=.5]
data_[1, results[1, :]>=.5] = 1.*(results[0, results[1, :]>=.5]==1)
data.append(data_)
for ax in [ax1, ax2]:
ax.axis([0., 1., -.1, 1.1])
_ = ax.set_xlabel('contrast')
_ = ax1.set_ylabel('detected lower')
_ = ax1.legend(['lower', 'upper'], loc='right')
_ = ax2.set_ylabel('detected upper')
_ = ax2.legend(['lower', 'upper'], loc='right')
(note the subplots are equivalent up to a flip)
One could not fit psychometric curves for this 2AFC. See for instance https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3110332/
In [30]:
n_hist = 10
C_bins = np.linspace(0.5, 1., n_hist)
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(111)
for data_ in data:
hist_correct, bins = np.histogram(data_[0, data_[1, : ]==1], bins=C_bins, normed=False)
hist_total, bins = np.histogram(data_[0, : ], bins=C_bins, normed=False)
plt.bar(bins[:-1] + (bins[1]-bins[0])/2, (1.*hist_correct)/(1.*hist_total), width=.95*(bins[1]-bins[0]), alpha=.1)
_ = ax.axis([0.5, 1., 0., 1.1])
_ = ax.set_xlabel('contrast')
let's explore another way:
In [31]:
import pypsignifit as psi
data : an array or a list of lists containing stimulus intensities in the first column, number of correct responses (nAFC) or number of YES- responses in the second column, and number of trials in the third column. Each row should correspond to one experimental block. In addition, the sequence of the rows is taken as the sequence of data aquisition. Alternatively, the relative frequencies of correct responses resp YES responses can be given.
In [32]:
stimulus_intensities, percent_correct = np.array([]), np.array([])
for data_ in data:
#print data_.shape
#stimulus_intensities = np.hstack((stimulus_intensities, data_[0, : ]))
#percent_correct = np.hstack((percent_correct, data_[1, : ]))
#print data_[0, : ].min(), data_[0, : ].max(), data_[1, : ].mean()
hist_correct, bins = np.histogram(data_[0, data_[1, : ]==1], bins=C_bins, normed=False)
hist_total, bins = np.histogram(data_[0, : ], bins=C_bins, normed=False)
stimulus_intensities = np.hstack((stimulus_intensities, bins[:-1] + (bins[1]-bins[0])/2))
percent_correct = np.hstack((percent_correct, (1.*hist_correct)/(1.*hist_total)))
num_observations = len(stimulus_intensities)
print 'num_observations: ', num_observations
psidata = np.vstack((stimulus_intensities, percent_correct, np.array([2]*num_observations)))
BootstrapInference
In [34]:
# see http://psignifit.sourceforge.net/DESIGN.html
nafc = 2
constraints = ( 'unconstrained', 'unconstrained', 'Beta(2,20)')
B = psi.BootstrapInference ( psidata, priors=constraints, nafc=nafc )
# (threshold, slope, lapse rate)
print B.estimate
How well do these parameters describe the data? The deviance is a measure that describes the goodness of fit for a model, based on the sum of the squares error metric:
In [35]:
print B.deviance, B.getThres()#, B.getCI(1)
In [36]:
_ = psi.psigniplot.plotPMF(B, xlabel_text='contrast', showaxes=True)
In [37]:
B.sample()
In [38]:
psi.GoodnessOfFit(B)
- Top left: displays the given data points and the fitted psychometric function. Thresholds and confidence intervals are plotted at three levels (default: 25%, 50% and 75% ). Mind that the y-axis starts at 0.5 (the guess rate in a 2AFC task), therefore the 50% threshold is located at . signifies the deviance value.
- Bottom left: histogram of bootstrapped deviance values (default = 2000 samples). The 95% confidence limits are indicated by red dotted lines and the actually observed deviance is indicated by the solid red line. If the observed deviance is outside the 95% confidence limits that indicates a bad fit and you will receive a warning message.
- Top middle: deviance residuals are plotted as a function of the predicted correct response rate of the model (x-axis corresponds to y-axis in panel 1). This plot helps you to detect systematic deviations between the model and the data. The dotted line is the best linear fit that relates deviance residuals to the predicted correct response rate. Rpd gives the numerical value of that correlation. Note that the residuals are scaled to account for differences in the variability of a binomially distributed random variable (e.g. maximum variance at p=0.5).
- Bottom middle: histogram of bootstrapped correlation coefficients for the correlation between residuals and performance level (same logic applies as in panel 2). Dotted lines denote 95% intervals of the sampled correlation coefficients, the solid line marks the observed correlation between model prediction and deviance residuals.
- Top right: deviance residuals are plotted as a function of block index i.e. the sequence in which the data were acquired (WARNING: this graph can be properly interpreted only when stimulus intensities were fixed in separate blocks). If the observer was learning, the fitted linear correlation between residuals and block index should be positive.
- Bottom right: histogram of bootstrapped correlation coefficients for the correlation between deviance residuals and block index (same logic applies as in panel 2 and 4).
In [13]:
priors = ( 'Gauss(0,5)', 'Gamma(1,3)', 'Beta(2,20)' )
In [14]:
mcmc = psi.BayesInference ( psidata, priors=priors, nafc=nafc )
In [15]:
mcmc.estimate
Out[15]:
In [16]:
psi.ConvergenceMCMC ( mcmc )
In [17]:
psi.GoodnessOfFit ( mcmc )
In [18]:
psi.ParameterPlot ( mcmc )
Out[18]:
In [19]:
psi.plotInfluential(mcmc)
In [19]:
In []:
Testing basic functions of Motion Clouds
Motion Clouds: raw principles¶
Motion Clouds are synthesized textures which aim at having similar characteristics as natural images but with controlled parameters. There are many ways to achieve these results and this notebook aims at showing that different procedures from different communities (neurioscience, modelling, computer vision, ...) may produce similar results.
In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
Using Fourier ("official Motion Clouds")¶
In [2]:
import sys
sys.path.append('..')
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
Using mixtures of images¶
In [3]:
import scipy.misc
lena = scipy.misc.lena() * 1.
lena -= lena.mean()
lena /= lena.std()
print lena.shape
In [4]:
plt.imshow(lena, cmap=plt.cm.gray)
Out[4]:
In [5]:
def noise(image=lena):
for axis in [0, 1]:
image = np.roll(image, np.random.randint(image.shape[axis]), axis=axis)
return image
In [6]:
plt.imshow(noise(), cmap=plt.cm.gray)
Out[6]:
In [7]:
plt.imshow(noise(), cmap=plt.cm.gray)
Out[7]:
Using ARMA processes¶
Now, we define the ARMA process as an averaging process with a certain time constant $\tau=30.$ (in frames).
In [8]:
def ARMA(image, tau=30.):
image = (1 - 1/tau)* image + 1/tau * noise()
return image
initializing
In [9]:
image = ARMA(lena)
plt.imshow(image, cmap=plt.cm.gray)
Out[9]:
In [10]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
Out[10]:
In [11]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
Out[11]:
In [12]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
Out[12]:
In [16]:
%cd..
N_frame = 1024
z = np.zeros((lena.shape[0], lena.shape[1], N_frame))
z[:, :, 0] = image
for i_frame in range(1, N_frame):
z[:, :, i_frame] = ARMA(z[:, :, i_frame-1])
mc.anim_save(.5 + .5*z, filename='results/arma')
In [18]:
mc.notebook = True
mc.in_show_video(name='arma')
Testing discrimination between Motion Clouds with psychopy
Analysis of a discrimination experiment¶
In the psychopy_discrimination.py script, we implemented an experiment to test whether the overall shape of motion clouds could influence discrimination of speed.
Herein, we analyse the data that was collected over different sessions and try to draw some conclusions.
In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
The experiment¶
In [2]:
#%pycat psychopy_discrimination.py
Analysis - version 1¶
In a first version of the experiment, we only stored the results in a numpy file.
In [11]:
!ls ../data/discriminating_*
In [12]:
!ls ../data/
In [4]:
import glob
for fn in glob.glob('../data/discriminating_*npy'):
results = np.load(fn)
print fn, results.shape
print('analyzing results')
print '# of trials :', np.abs(results).sum(axis=1)
print 'average results: ', (results.sum(axis=1)/np.abs(results).sum(axis=1)*.5+.5)
In [5]:
%whos
Another solution is to use:
In [6]:
#data= 1
#np.savez('toto', data=data, results=results)
#print np.load('toto.npz')['data']
or directly psychopy.misc:
In [7]:
from psychopy import misc
In [9]:
#info = misc.fromFile('../data/discriminating.pickle')
In [10]:
alphas = [-1., -.5, 0., 0.5, 1., 1.5, 2.]
fileName = '../data/discriminating_' + info['observer'] + '.pickle'
info['alphas'] = alphas
info['results'] = results
#numpy.savez(fileName, results=results, alphas=alphas)
misc.toFile(fileName, info)
Analysis - version 2¶
In the second version, we stored more data:
In [14]:
import glob
from psychopy import misc
fig = plt.figure()
ax = fig.add_subplot(111)
for fn in glob.glob('../data/discriminating_v2_*pickle'):
data = misc.fromFile(fn)
print fn, results.shape
print('analyzing results')
alphas = np.array(data['alphas'])
print ' alphas = ', alphas
print '# of trials :', np.abs(data['results']).sum(axis=1)
print 'average results: ', (data['results'].sum(axis=1)/np.abs(data['results']).sum(axis=1)*.5+.5)
width = (alphas.max()-alphas.min())/len(alphas)
ax.bar(data['alphas'] - width/2, (data['results'].sum(axis=1)/np.abs(data['results']).sum(axis=1)*.5+.5), width=width, alpha=.3 )
# for i_trial in xrange(data['results'].shape[1]):
# ax.scatter(data['alphas'][data['results']
ax.set_xlabel('alpha')
ax.set_ylabel('% correct')
Out[14]:
In []:
# TODO : correct for the number of trials / make global reading / make a new draw of alphas
# TODO : test f_0
In []: